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)  The  development  of  a  consistent  penalty  -  Galerkin  finite  element  algorithm 
for  small-scale  atmospheric  flows  over  complex  terrain  is  summarized.  Various 
aspects  of  the  development  such  as  pressure  smoothing,  the  development  of  an 
appropriate  variable  step  time  Integrator,  and  appropriate  outflow  boundary 
conditions  are  discussed  and  pertinent  references  cited.  Illustrative 
simulations  are  presented. 
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PROBLEH  DESCRIPTION 

Introduction 

There  are  many  problems  of  interest  to  fluid  dynamicists  as  well  as 
meteorologists  which  require  a  small-scale,  or  limited  area,  model  of  flow  and 
possibly  concomitant  transport  of  energy  and  species  over  irregular  terrain. 
For  example  stratified  flow  over  hills,  or  cities,  and  the  associated  transport  of 
pollutants  or  atmospheric  boundary  layer  studies  in  which  the  resolvable  hor¬ 
izontal  and  vertical  scales  are  of  the  same  order  of  magnitude  and  the  hydros¬ 
tatic  approximation  is  invalid.  Considerable  progress  has  been  made  in 
developing  numerical  techniques  to  simulate  models  of  stratified  flows  without 
topography  but  the  inclusion  of  general  topographic  effects  has  only  recently 
received  a  thorough  consideration  from  a  numerical  simulation,  as  well  as 
operational  standpoint,  by,  for  example,  Gal-Chen  and  Somerville  (1975),  Clark 
(1977),  Mason  and  Sykes  (1979  a,  b).  Mahrer  and  Pielke  (1975,  1977)  and  Hauss- 
ling  (1977).  The  latter  authors  employ  a  terrain-following  coordinate  system 
which  maps  the  irregular  domain  into  a  rectangular  domain  and  then  develops 
finite  difference  stencils  for  the  various  terms  appearing  in  the  transformed 
balance  equations  which  must  be  satisfied  by  the  primitive  variables.  Many  of 
the  aspects  of  these  modeling  efforts  are  reviewed  in  Pielke  (1981). 

The  main  thrust  of  our  ARO  sponsored  research  effort  is  the  development 
and  assessment  of  a  primitive  variable,  non-hydrostatic  Galerkin-finite  element 
technique  which  can  accommodate  general  topography,  is  capable  of  automati¬ 
cally  generating  a  spatial  discretization  which  conserves  certain  global  quanti¬ 
ties  if  desired  (Lee,  Gresho,  Chan,  Sani  and  Cullen  (1981)  )  and  is  robust  enough 
to  handle  very  complex  flows.  We  have  developed  and  assessed  the  potential  of 
a  Galerkin  finite  element  penalty  technique  for  flow  and  transport  problems, 
and,  I  believe,  have  clearly  demonstrated  its  potential  in  satisfying  the  above 
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r  squire  me  nis. 

BASIC  DESCRIPTION  OF  PROTOTYPE  SYSTEM  AND  TECHNIQUE 

In  order  to  simplify  the  discussion  of  the  important- features  of  the  pro¬ 
posed  technique  the  anelastic  model  used  by  Gal-Chan  and  Somerville  (1975) 
will  be  used.  (More  complicated  anelastic  models  with  higher  order  turbulence 
closure  schemes  are  compatible  with  the  technique  but  of  questionable  utility 
in  3-D  simulation  on  current  generation  computers.)  The  basic  mathematical 
model  is: 

V  (p0u)  =  0  (1) 

0(pou) 

— gj +  V  (p0uu)  =  -  Vp  +  pg  +  V  T  (2) 

d(p08>) 

+  V-Gw)  =  VH  (3) 

Here  t  and  H  are  respectively,  the  stress  tensor  and  heat  flux  which,  in  general, 
are  nonlinear  functions  of  the  local  state  of  the  flow.  The  Coriolis  effect  has 
been  neglected  but  is  easily  included. 

In  applying  a  finite  difference  (FD)  method  these  equations  are  usually 
transformed  to  a  terrain  following  coordinate  system  and  a  discretized  version 
is  generated.  In  effecting  the  latter  oftentimes  a  diagnostic  pressure  equation 
is  generated  and  also  discretized.  In  many  cases  the  discretized  version  of  the 
diagnostic  pressure  equation  is  rather  complicated  and  its  solution  involves 
some  sort  of  iteration  scheme.  One  of  the  interesting  features  of  the  penalty 
Galer kin- Finite  element  (PGFEM)  technique  proposed  here  is  that  the  basic 
discretized  system  is  one  which  only  involves  nodal  velocities  and  the  nodal 
pressures  are  obtained  by  post-processing  of  the  nodal  velocities.  The  basic 
idea  appears  to  be  due  to  Lions  (1972)  and  belongs  to  the  same  family  of 
slightly  compressible  approximations  first  introduced  by  Chorin  (1967)  but  is 
not  identical  and  is  implemented  here  in  a  completely  different  manner.  The 


standard  GFEM  approach  to  anelastic  fluid  mechanical  systems  would  involve 
interpolating  both  the  velocity  and  pressure  fields  and  generating  nodal  equa¬ 
tions  as  outlined  in  a  subsequent  section.  In  a  mixed  scheme  both  nodal  pres¬ 
sures  and  velocities  appear  as  unknowns,  not  necessarily  at  each  node  since  it 
is  necessary  to  interpolate  the  pressure  field  with  a  lower  order  polynomial 
than  the  velocity  field  to  obtain  a  well-posed  algebraic  problem  for  the  nodal 
unknowns.  (See.  for  example,  Taylor  and  Hood.  1973;  and  Sani.  Gresho  and  Lee 
(1980),  Sani.  Gresho,  Lee  and  Griffiths  (1981a),  Sani,  Gresho.  Lee.  Griffiths  and 
Engelman  (1981b). 

In  order  to  introduce  the  proposed  penalty  GFE.  i.e.  PGFE.  technique  one 
introduces  the  following  perturbed  system. 


7(poU«)  =  -ept, 

0(Po««)  >  1 

oi —  ♦  -  jjV  (poU,)|u« 


(4) 


*  -  yp«  +  pg  +  V  Tt ,  (5) 

8(Po«0  .  i 

—ft —  +  V-GxjIMp,)  -  (poU,)^,  =  V  H,  (6) 

in  which  e,  the  "penalty  parameter",  appears  as  a  small  perturbation  parame¬ 
ter.  By  using  the  results  of  Teman  (1988,  1976)  and  placing  some  modest  res¬ 
trictions  on  the  form  of  t  and  H  it  appears  possible  to  formally  establish  that 
the  two  systems  differ  by  an  0  (t)  error  which,  of  course,  implies  that  the  two 
systems  (l)-(3)  and  (4)-(8),  are  identical  as  e-*0.  The  addition  of  terms  involving 

g-7  (p0u«)  insure  that  the  discretized  version  of  the  advection  terms  is  conser- 

‘▼ative  (up  to  time  truncation  errors)  and  lions  (1972)  has  shown  that  these 
terms  can  be  added  without  jeopardizing  the  accuracy  of  the  method  as  long  as 
e  is  sufficiently  small.  The  implication  of  these  terms  has  been  assessed  in  one 
aspect  of  the  work  initiated  under  the  current  research  effort  Lee,  Gresho, 
Chan,  Sani  and  Cullen  (1981)).  It  was  found  that  the  GFE,  or  the  PGFE,  method 
is  relatively  insensitive  (compared  to  a  FD  method)  to  the  inclusion  of  these 
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terms  to  obtain  a  flux  divergence  form  of  the  equations;  however,  the  computa¬ 
tion  of  derived  quantities  can  be  greatly  effected  (Gresho,  Lee  and  Sani  (1981)) 
and  hence  we  routinely  employ  the  flux  divergence  form. 

The  PGFE  method  utilized  in  our  current  algorithm  leads  to  a  formulation 
in  which  the  nodal  pressure  variables  do  not  explicitly  appear  —  a  cost  effective 
version  of  the  general  technique.  It  is  noteworthy  that  our  implementation 
does  not  formally  replace  p  before  discretizing  via  the  GFE  which,  as  is  well- 

known,  generates  a  "penalty  term",  — V(V  u)  which  must  be  "underintegrated” 

£ 

when  discretized  by  Galerkin's  method,  but  directly  discretizes  the  perturbed 
continuum  equations  via  GFE  using  a  C~x  finite  element  pressure  representa¬ 
tion.  The  latter  results  in  a  discretized  approximation  of  the  perturbed  con¬ 
tinuity  equation  which  can  be  solved  for  nodal  pressures  at  element  level,  a 
very  efficient  process.  Whenever  the  nodal  pressures  within  an  element  are 
desired  they  can  be  obtained  by  appropriate  post-processing  of  the  nodal  velo¬ 
cities  associated  with  the  element.  This  PGFE  method  is  referred  to  as  a  "con¬ 
sistent"  penalty  technique  and  its  implementation  and  advantages  are  dis¬ 
cussed  in  Engelman,  Sani,  Gresho  and  Bercovier  (1981). 

M 

In  the  implementation  of  the  scheme  the  approximate  solution  (u.p)  is 
represented  in  the  form 

u(*.0»2>/(0p/(*).  p(x,t)~'£pi(t)i//j(x)  (7) 

/■i  }• i 

The  finite  element  basis  set  {ty(x)|  (jV'j(x)P  is  nonzero  only  locally  over  a  few 
elements  (an  element),  i.e.,  is  of  local  support,  and  ip}  (V'>)  is  unity  at  node  j  and 
zero  at  other  nodes.  (In  certain  cases  it  is  convenient  to  relax  the  latter  pro¬ 
perty  for  fj.)  The  basis  $ty(x){  are  required  to  be  at  least  Ca  and  the  pressure 
basis  H'j(x)]  is  required  to  be  C~l  for  efficient  implementation  of  the  tech¬ 
nique.  Both  these  choices  are  compatible  with  the  weak  form  of  the  equations 
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associated  with  our  Galerkin  formulation  which  leads  to  a  set  of  coupled  first 
order  ordinary  differential  equation  of  the  form 

H^-  +  (  ~)Bu  +  Ku  +  A(u)u  =  b  (8) 

in  which  the  matrices  are  sparse  and  usually  banded  but,  in  general,  the 
bandwidths  are  much  larlarger  than  the  corresponding  system  generated  by  a 
regular  finite  difference  stencil.  This  property  implies  that  for  the  PGFEM  to  be 
competitive  with  a  FD  method  fewer  nodal  unknowns  must  be  used  to  achieve  a 
specified  accuracy  and  efficient  storage  and  solution  techniques  must  be  util¬ 
ized.  The  current  state  of  GFEM  (or  PGFEM)  development  is  such  that  very  few 
documented  comparisons  between  the  GFEM  and  FD  schemes  Eire  available  in 
the  literature  but  some  preliminary  coarse  mesh  comparisons  or  the  schemes  in 
simulating  flow  and  advection- diffusion  were  encouraging  (Gresho,  Lee  and  Sani 
(1980)),  The  cost-effectiveness  is  important  and  an  area  of  current  research; 
however,  robustness  is  also  a  very  important  attribute  and  a  recent  test  prob¬ 
lem  (de  Vahl  Davis  and  Jones  (1981))  points  out  the  superiority  of  Galerkin  fin¬ 
ite  element  techniques  in  this  regard. 

The  appearance  of  the  non-diagonal  mass  matrix  II  in  (8)  implies  the 
implicit  character  of  the  technique;  in  non-penalty  Galerkin  finite  element 
techniques  a  so-called  lumped  version  which  replaces  H  by  a  diagonal  matrix 

can  often  be  effectively  used  but  here  with  (—  )>>1  an  explicit  scheme  is  inap- 

c 

propriate  because  of  very  stringent  stability  limitations,  Le.,  Af~0(e)  and  typi¬ 
cally  e«-10"e  -  10-8.  In  fact,  there  is  a  false  transient  induced  by  the  method  (a 
ramification  of  the  singularity  of  B  and  smallness  of  e)  which  must  be  dealt  with 
properly.  These  points  are  discussed  in  more  detail  in  the  Summary  of  Results 
Section. 
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SUHHAKY  OF  RESULTS 

The  emphasis  in  the  program  from  its  inception  has  been  one  of  develop¬ 
ment  coordinated  with  theory.  Our  initial  study  focused  firstly,  on  understand¬ 
ing  and  implementing  the  penalty  technique  and  secondly,  on  the  ability  of  the 
technique  to  generate  accurate  solutions  to  difficult  problems.  During  the 
course  of  the  study  detailed  investigations  were  carried  out  on  (l)  spurious 
pressures  and  consistent  penalty  in  simulation  of  incompressible  media;  (2)  the 
"false  transient"  associated  with  the  consistent  penalty  technique  and  the 
design  of  a  variable  step  time  integration  with  local  error  control  to  cope  with 
this  problem  in  a  cost-effective  manner;  (3)  the  generation  and  evaluation  of 
new  elements;  (4)  the  derivation  of  "consistently"  derived  normals  on  the  boun¬ 
dary  of  a  computational  domain  which  is  represented  in  the  GFE  method  by  a 
C°  representation;  (5)  initial  studies  of  turbulence  parameterization;  (6)  the 
development  of  pre-  and  post-  processing  techniques;  (7)  benchmark  simula¬ 
tions  assessing  and  illustrating  features  of  the  technique.  Since  these  studies 
have  formed  the  basis  of  the  published  manuscripts  cited  in  the  Publications 
Section  and  are  readily  available,  only  the  highlights  of  certain  aspects  are 
presented  here. 

1.  Spurious  Pressures  and  Consistent  Penalty  formulation 

Discretized  approximations  to  the  incompressible  Navier-Stokes  equations, 
in  the  primitive  variable  (velocity-pressure)  formulation,  especially  when  gen¬ 
erated  via  the  Galerkin  finite  element  method  (GFEM),  have  been  plagued  with 
confusion  regarding  the  ’appropriate’  workable  combinations  of  velocity  and 
pressure  approximations.  Early  on  it  was  discovered  by  Hood  and  Taylor  (1973) 
that  equal-order  Interpolation  on  conforming  quadrilateral  elements,  wherein 
the  same  basis  functions  are  used  for  representing  velocity  and  pressure, 
causes  difficulty  in  the  pressure  solution.  They  obtained  better  results  when 
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using  mixed  interpolation  (the  basis  functions  for  pressure  were  one  order 
lower  than  those  for  velocity),  and  suggested  an  explanation  in  the  form  of 
.  "balancing  residuals'*  from  momentum  and  continuity  equations.  Their  expla¬ 
nation,  although  intuitively  appealing,  was  judged  inadequate  by  Olson  and 
Tuann  (1978)  who  explained  the  results  in  terms  of  the  eigenvalues  of  a  single 
element  (equal  interpolation  always  generated  one  or  more  eigenvectors  which 
contained  only  pressure  and  corresponded  to  zero  eigenvalues;  they  claimed 
that  these  were  spurious  and  were  the  cause  of  the  failure). 

Even  when  mixed  interpolation  is  employed,  however,  there  are  cases 
where  numerical  difficulties  are  encountered.  In  particular,  the  simplest  possi¬ 
ble  approximation  employs  piecewise  linear  approximation  for  velocity  (bilinear 
on  quadrilaterals)  and  piecewise  constant  approximation  for  pressure.  This  ele¬ 
ment  has  been  found  to  work  well  in  some  cases  and  poorly  in  others;  in  solid 
mechanics,  see  Argyris  et  al  (1974)  and  Nagtegaal,  Parks  and  Rice  (1974),  and  in 
fluid  mechanics,  see  Fabayo  (1977),  Huyakorn  et  al  (1978),  Hughes,  Liu  and 
Brooks,  (1979),  and  Lee.  Gresho  and  Sani  (1979).  For  certain  combinations  of 
boundary  conditions  and  element  distributions  over  a  domain,  the  solutions 
display  acceptable  velocities  but  totally  spurious  pressures;  the  pressures  are 
afflicted  with  the  "checkerboard  syndrome,"  wherein  they  display  oscillations 
which  are  frequently  of  one  sign  on  all  ’black’  squares  and  of  the  opposite  sign 
on  all  ‘red*  squares.  These  pressure  patterns  have  also  been  encountered  using 
certain  finite  difference  discretization  techniques  (Fortin  (1972),  Pracht  and 
Blackbill  (1978)  so  that  the  affliction  is  not  intrinsic  to  GFEM  formulations;  for 
-  example,  Chorin  (1969)  has  encountered  four  spurious  pressure  patterns  when 
solving  a  consistently  centered  difference  approximation  to  the  Navier-Stokes 
equations  in  two-dimensions  (eight  in  3-D).  Finally,  similar  behavior  can  also 
occur  using  higher-order  elements;  e.g.  the  quadrilateral  element  with  biqua- 


T 
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dratic  velocity  (9-node)  and  bilinear  pressure  with  nodes  at  the  2x2  Gauss 
points  can  display  a  checkerboard  syndrome. 

In  the  related  penalty  method  (  PGFEM  )  which  is  of  interest  herein,  the 
solenoidal  constraint  on  the  velocity  field  is  slightly  weakened,  to  p  =  -  \7  u, 

where  X  =  (  — )  is  the  penalty  parameter.  The  results  from  the  penalty  approxi- 
c 

mation  will  be  close  to  those  using  (lb)  if  the  following  inequalities  are  satisfied: 
e<  <1,  where  e  is  the  "unit  roundoff  level”  of  the  computer:  for  example, 

on  our  CDC-720,  e«10~u  and  the  penalty  method  'works  well’  for 

«108  <  X/pi  <  «10®.  With  perfect  arithmetic  (s  =  0),  the  penalty  results  would 
converge  to  those  from  (l)  as  X— >■».  For  large  but  finite  X,  the  results  should 
typically  differ  by  0(1/ X):  actually,  we  usually  observe  «6  digits  of  agreement 
for  10®  ,<  X  10®.  Since  the  pressures  obtained  from  the  penalty  method  can 
also  be  plagued  with  the  checkerboard  syndrome,  it  is  appropriate  to  consider 
it  also. 

It  is  shown  in  Sani  et  al  (1981  a,  b)  that  the  spurious  pressures  are  ramifi¬ 
cations  of  non-physical  elements  in  the  null  space  of  the  discretized  pradient 
operator,  i.e.,  solution  to  the  GFEM  discretized  form  of  (4)  -  (5)  is  of  the  form 
(u.p)  =  (0 ,p)  where  p  is  not  a  constant  vector.  The  consequences  and  filtering 
of  these  spurious  pressures  on  regular  as  well  as  isoparametric  meshes  is  dis¬ 
cussed  in  the  above  references.  As  alluded  to  above  the  penalty  method  advo¬ 
cated  herein  is  theoretically  devoid  of  spurious  pressures  but,  in  fact,  finite 
precision  numerical  simulations  using  certain  elements  can  exhibit  oscillations 
but  usually  of  very  small  amplitude.  If  desired  these  oscillations  may  also  be  fil¬ 
tered  by  techniques  alluded  to  above.  However,  in  most  cases  one  can  select 
velocity-pressure  approximations  which  do  not  exhibit  such  spurious  pressures, 
and  for  example,  in  2D  the  biquadratic  velocity-linear  pressure  interpolant 
which  we  reported  in  Engelman  et  al  (1982)  appears  to  be  an  optimal 


combination. 


An  example  comparing  "filtered  pressures"  from  a  biquadratic  C°  velocity  - 
linear  C-'  pressure  element  mixed  mode  simulation  with  the  unfiltered  pres¬ 
sures  from  the  consistent  penalty  approach  advocated  herein  using  the  same 
element  pair  as  well  as  a  bilinear  C°  velocity  -  piecewise  constant  C-1  pressure 
element  is  illustrated  in  Figure  1. 


A  Driven  Cavity  Experiment 


Since  the  mesh  is  symmetric,  and  the  "penalty"  pressures  satisfy  f  pK  =  0 

n 

only  one-half  of  the  pressures  are  displayed;  the  pressures  at  the  Gauss  points 
in  each  element  shown  on  the  left  side  are  the  negative  of  those  printed  on  the 
right  side.  For  the  mixed  interpolation  case,  the  spurious  pressure  was  first  fii- 

a 

* 

tered,  then  the  hydrostatic  level  was  adjusted  to  agree  with  that  from  the 
penalty  case.  The  first  entry  at  each  "node"  is  that  from  mixed  interpolation; 
the  second  entry  is  the  9-node  penalty  pressure,  and  the  bottom  entry  is  the 
4-node  penalty  pressure. 

It  is  apparent  in  Figure  1  that  all  pressures  look  fairly  reasonable  in  the 
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high  pressure  region  and,  considering  the  coarseness  of  the  mesh,  agree  quite 
well  with  each  other.  The  differences  show  up  in  the  low  pressure  region 
wherein  a  ''lingering  local  pressure  oscillation  is  present  in  both  penalty 
results.  Finally,  as  mentioned  earlier,  if  the  9-node  penalty  pressures  are  fil¬ 
tered,  the  resulting  pressures  agree  with  the  filtered  mixed  interpolation  pres¬ 
sures  to  *0(1/ X)  and  are  the  physical  pressures  in  the  sense  that  they  are  the 
only  ones  "seen"  in  the  momentum  equations.  The  latter  pressure  would  have 
been  obtained  without  filtering  by  employing  a  linear  C~x  pressure  interpolant 
initially  -  a  ramification  of  the  optimality  of  this  element. 

A  detailed  theoretical  discussion  of  the  filtering  schemes  as  well  as  addi¬ 
tional  numerical  examples  are  presented  in  Lee,  Gresho  and  Sani  (1979);  more- 
overs  a  detailed  discussion  of  the  Galerkin-finite  element  penalty  method  is 
presented.  It  is  noteworthy  that  the  implementation  of  the  technique  herein  is 
via  the  "consistent  penalty  method"  and  not  the  popular  under-integration 
implementation  most  often  employed  in  the  engineering  literature.  While  the 
two  techniques  may  be  equivalent  on  a  rectangular  mesh  there  can  be  signifi¬ 
cant  differences  on  isoparametric  meshes  composed  of  generalized  quadrila¬ 
terals.  The  latter  is  illustrated  in  Table  1,  detailed  in  Engelman,  Sani,  Gresho 
and  Bercovier  (1982)  and  represented  an  important  contribution  of  our  initial 
investigation.  One  of  the  benchmark  simulations  used  in  generating  a  com¬ 
parison  of  the  methods  was  a  "colliding  flow".  In  this  flow,  fluids  flow  from  +y„ 
and  — ym  to  collide  at  y  =  0.  The  solution  over  the  entire  x-y  plane  is 

u  t  -  Cx 
uy  =  -Cy 

p  =  -fp(*8  +  ye) 

The  region  Oirsl.O^y^l  was  modeled  with  symmetry  planes  (zero  normal  velo¬ 
city  and  tangential  stress)  prescribed  along  x  -  0  and  y  =  0  and  the  exact  solu- 
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tion  along  x  -  1  and  y  =  1.  The  parameters  used  were  C  =  4,  density  p  =  1  and 
viscosity  p,  =  .0667. 

The  approximate  solution  to  each  problem  was  obtained  using  both  the 
reduced  integration  and  the  consistent  construction  of  the  penalty  matrix  for 
ptQ j  and  piP\  on  three  different  meshes.  The  meshes  were:  a  regular  6x6 
mesh  of  nine  node  squares,  a  6x6  mesh  of  irregular  quadrilate rials,  and  a  8x6 
mesh  of  curved  isoparametric  elements.  The  isoparametric  mesh  was  generated 
by  moving  each  mid-side  node  inwards  or  outwards  a  distance  of  5%  of  the  ele¬ 
ment  width.  The  results  of  the  simulation  are  presented  in  Table  1.  They  are  in 
complete  agreement  with  the  theoretical  predictions;  that  is,  on  both  the  regu¬ 
lar  and  distorted  quadrilateral  meshes  both  the  reduced  and  consistent  con¬ 
structions  of  the  penalty  matrices  give  essentially  identical  solution  errors.  It 
is  only  on  the  curved  isoparametric  meshes  that  the  differences  between  the 
reduced  and  consistent  formulations  manifest  themselves.  Here  we  see  a  large 
error  in  the  pressure  solution  when  reduced  integration  is  used  which 
decreases  markedly  when  a  consistent  construction  is  employed.  Only  the 
results  for  the  isoparametric  mesh  are  reported  in  Table  1.  as  the  regular  and 
quadrilateral  meshes  gave  nearly  identical  errors  for  all  three  formulations. 
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Mesh 

Pressure 

Formulation 

Rel.  Error 
Velocity  • 

Rel.  Error 
Pressure 

Isoparametric 

Elements 

Reduced  Integration 
P«9i 

.05057 

.02384 

Consistent  p  e  Q  j 

.03185 

.02129 

Consistent;) 

.01876 

.01601 

TabLe  1:  2-D  Colliding  Flow  Problem 

2.  False  Transient  Associated  With  The  Penalty  Technique. 

The  penalty  method  used  in  our  small  scale  atmospheric  model  weakens 
the  solenoidal  constraint  (l)  to  the  weaker  one  (4).  As  a  consequence,  the 
infinite  speed  of  propagation  of  a  pressure  signal  in  an  incompressible  fluid  is 
replaced  by  an  0(A)  "diffusive"  speed  of  propagation  -  a  property  which  can  lead 
to  a  non-physical  transient.  Here  A=l/s  and  e<  <1.  Since  the  penalty  parame¬ 
ter  is  large,  the  initial  nonphysical  transient  is  fast,  occurring  on  a  time  scale  of 
0(1/ A),  and  reflects  the  slight  compressibility  inherent  in  the  weakened 
solenoidal  constraint.  While  the  analogy  with  a  "special”  slightly  incompressible 
fluid  has  been  noted  by  others,  it  appears  that  the  associated  non-physical 
transient  has  received  little  attention  in  the  literature.  In  Sani  et  al  (1981c) 
attention  focused  on  the  existence,  characterization  and  remedy  of  the  initial 
-  non-physical  portion  of  the  transient  associated  with  a  well-posed  (in  the 
incompressible  sense)  problem  and  also  illustrate  the  physically  uninteresting 
transient  associated  with  an  ill-posed  problem.  The  nature  of  the  false  tran¬ 
sient  can  be  demonstrated  by  considering  the  development  of  a  Poisieulle  flow. 
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Since  the  penalty  parameter,  X,  is  very  large,  it  can  be  shown  that  an  expli¬ 
cit  time  integration  algorithm  leads  to  time  steps  (Af)  of  0(1/ X)  for  stability  - 
an  unacceptable  restriction  imposed  by  the  non-physical  part  of  the  transient. 
Even  implicit  algorithms  which  are  non- dissipative  can  be  severely  limited  in 
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time  step  size;  for  example,  if  the  trapezoid  rule  algorithm  is  chosen,  undamped 
oscillations  will  occur  if  XAf  is  much  greater  than  unity.  Consequently,  a  time 
integrator  was  designed  with  the  capability  of  starting  out  in  a  constant,  or 
variable,  Af  backward  Euler  (BE)  mode  (an  implicit,  first-order,  dissipative 
scheme)  and  switching  to  a  variable  Af  BE  mode  or  trapezoid  (TR)  mode  (an 
implicit,  second-order,  non-dissipative  scheme)  after  a  preassigned  number  of 
time  steps.  Using  this  algorithm  one  can  either  follow  the  entire  transient  to 
within  a  prescribed  local  time  truncation  accuracy  via  a  variable  step  method 
exclusively,  or  follow  only  the  physical  part  of  the  transient  by  initially  taking  a 
few  BE  steps  at  a  constant  Af»(l/X)  and  then  switching  to  variable  step  BE  or 
TR  for  efficiency.  The  initial  non-physical  portion  of  the  transient  appears  as 
"high  frequency  signals"  to  the  BE  scheme  and  is  therefore  quickly  damped  by 
the  numerical  diffusion  associated  with  this  scheme.  Consequently,  this  stable 
time  integration  strategy  can  over-look  the  initial  non-physical  transient  while 
using  reasonable-sized  time  steps,  i.e.,  time  steps  not  0(c),  and  can  be  up  to 
second  order  accurate  in  the  remaining  (physical)  portion  of  the  transient 
while  automatically  varying  Af,  based  solely  on  temporal  accuracy  require¬ 
ments.  The  details  are  reported  in  Sani  et  al  (1981c)  and  the  basic  ideas  can  be 
exploited  by  others  in  other  applications. 

3.  Consistently  Derived  Boundary  Normals 

It  is  often  required  to  specify  normal,  or  tangential,  velocities  along  the 
boundary  of  a  domain,  or  equivalently  along  the  element  boundaries  forming 
the  finite  element  covering  of  a  domain.  In  general,  these  boundaries  are  not 
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along  level  surfaces  of  the  coordinate  system  for  the  problem.  For  example,  in 
modeling  flows  with  free  surfaces  such  as  the  ocean-atmosphere  interaction  or 
in  raising  the  computational  domain  "off  the  ground"  when  assuming  the  lowest 
level  flow  is  well  represented  by  a  locally  equilibrium  log-linear  profile.  We 
addressed  this  issue  in  Engelman,  Sani  and  Gresho  (1982)  where  a  detailed 
analysis  is  presented  and  it  is  shown  that  if  the  normal  vector  computation  is 
not  done  in  a  manner  consistent  with  the  incompressibility  constraint,  errone¬ 
ous  results,  especially  in  pressure,  can  result. 

4.  Outflow  Boundary  Conditions 

During  our  study  it  was  demonstrated  mainly  by  numerical  experiment 
that  appropriate  computational  outflow  boundary  conditions  for  neutral  and 
special  stratified  flows  could  be  obtained  by  using  the  "natural"  boundary  con¬ 
ditions  associated  with  the  Galerkin  formulation,  i.e.,  constant  normal,  fn,  and 
in  some  cases  tangential  tt,  stress  on  an  outflow  boundary.  (Here  fn  =  nn:r  and 
f*  =  tar.)  While  these  conditions  are  not  generally  physically  correct  for  the 
continuum  system  (outflow  boundaries  imposed  by  the  requirements  of  a 
bounded  computational  domain  are  unphysical  and  apriori  the  "correct"  boun¬ 
dary  conditions  are  unknown)  benchmark  numerical  experiments  indicate  that 
such  boundary  conditions  lead  to  accurate  interior  numerical  results  especially 
in  advection  dominated  cases.  For  example,  in  the  simulation  of  the  shedding 
of  a  von  Karmen  vortex  street  at  Re  =  110  shown  in  Figure  4,  the  outflow  boun¬ 
dary  conditions  for  this  neutral  flow  were  fn  =  0,  f*  =  0.  Notice  that  the  vortex 

A 

0 

street  leaves  the  computational  domain  without  noticeable  mesh  effects. 
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Fig.  4.  Relative  Streamlines  As  Seen 
By  An  Observer  Moving  With 
The  Nominal  Fluid  Velocity 

It  is  also  noteworthy  that  the  weak  form  of  the  boundary  constraint  fn  =  const, 
is  not  enforced  locally  but  in  a  least  square  error  sense  over  the  outflow  boun- 
dary^only  in  the  limit  of  zero  mesh  size  is  the  constraint  imposed  locally.  It  is 
the  latter  feature  that  makes  this  natural  boundary  condition  which  allows  u 
and  p  to  be  computed  on  the  outflow  boundary  a  good  computational  outflow 
boundary  condition.  The  variation  of  fB  and  -p  on  the  outflow  boundary  of  the 
vortex  shedding  numerical  experiment  at  a  specific  time  is  displayed  in  Figure 
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Fig.  5.  fn  and  -p  On  Outflow  Boundary 
At  t  *  536.49 

It  is  noteworthy  that  fn«-p  since  on  this  outflow  boundary 

fru 

fn  =  -p  +2/x  — 

which  shows  that  even  at  Re  =  110  the  viscous  effects  at  the  outflow  boundary 
are  small  and  also  that  in  a  strong  sense  fn  *  const,  which  would  be  inappropri¬ 
ate  physically. 

While  the  solution  to  the  problem  of  appropriate  outflow  boundary  condi¬ 
tions  appears  to  be  reasonably  resolved  for  neutral  flows  with  small  Coriolis 
effects,  the  solution  in  the  case  of  many  stratified  flows  and/or  flows  with  signi¬ 
ficant  Coriolis  effects  requires  a  more  careful  analysis.  It  is  easily  demonstrated 
by  exact  solutions  as  well  as  numerical  experiments  that  setting,  for  example, 
fn  =  const.,  in  such  flows  leads  to  the  solution  of  a  different  problem  than  origi¬ 
nally  desired.  The  latter  is  a  ramification  of  both  velocity  and  pressure  field 
contributions  in  the  normal  stress  fn  coupled  with  the  occurrence  on  the  out¬ 
flow  boundary  of  a  significant,  apriori  unknown,  pressure  contribution  due 
mainly  to  hydrostatic  pressure  variations  in  a  stratified  flow  and/or  pressure 
variations  due  to  Coriolis  effects  in  a  rotating  flow.  (Centrifugal  can  be  handled 
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by  a  redefinition  of  the  pressure  variable  to  include  such  effects.) 


The  solution  to  this  outflow  boundary  condition  problem  is  being  dealt  with 
by  using  a  natural  boundary  of  the  form 

iii-r —  =  const, 
on 

which  again  is  not  physically  correct  pointwise  but  is  only  satisfied  in  a  weak 
sense  and  leads  to  an  acceptable  computational  outflow  boundary  condition 
except  in  the  case  of  a  piecewise  constant  pressure  interpolant  where  a  special 
fn  updating  procedure  must  be  used.  While  the  formal  implementation  of  the 
boundary  condition  is  straightforward  there  are  certain  theoretical  as  well  as 
computational  question  which  must  be  addressed  in  order  that  a  firm  basis  for 
its  use  is  established;  these  issues  are  addressed  in  Engelman,  Sani  and  Gresho 
(1982)  and  form  the  basis  of  some  ongoing  research.  Its  use  will  be  illustrated 
subsequently  in  one  of  the  examples  in  the  Benchmark  Simulation  Section. 

5.  Benchmark  Simulations 

During  the  course  of  this  study  a  number  of  numerical  experiments  were 
performed  in  order  to  check  the  basic  algorithm  as  well  as  to  judge  the  viability, 
effectiveness  and  robustness  of  new  techniques  during  the  course  of  their 
implementation  and  development.  These  simulations  varied  from  the  simplest 
for  which  the  algorithm  was  capable  of  generating  an  exact  solution  to  compli¬ 
cated  transient  stratified  flows  in  complex  domains.  The  results  of  some  of 
these  numerical  experiments  are  presented  here  in  order  to  illustrate  the 
capability  of  the  algorithm  as  well  as  some  of  the  pre-  and  post-processing 
techniques  which  were  either  developed  during  this  study  or  scavenged  from 
existing  algorithms. 

In  all  the  following  examples  illustrating  the  capabilities  of  the  algorithm  a 
C°  biquadratic  velocity  -  C~l  linear  pressure  element,  i.e.,  interpolation,  is 


employed  in  the  PGFEM  using  e  a  10"*.  Isoparametric  elements  were  used  to 
grade  the  mesh  and  to  track  the  complex  geometry.  (Each  element  appearing 
in  the  following  figures  illustrating  the  mesh  contain  18  velocity  degrees  of 
freedom  and  3  pressure  degrees  of  freedom.) 

1.  Steady,  Slow  Neutral  Flow  Over  A  Circular  Trench 


Fig.  6a.  Mesh 


Fig.  8c.  Streamlines 


Fig.  6d.  Isobars 

The  noteworthy  features  here  are  (i)  the  ability  to  easily  track  geometry 
and  grade  the  mesh  to  track  the  physics  if  known  apriori;  (ii)  the  quality  of  the 
solution  in  that  it  accurately  predicts  the  recirculation  zone  in  the  bottom  of 
the  trench,  a  feature  whose  existence  can  be  predicted  theoretically  and  (iii) 
the  fact  that  the  pressure  sigularities  which  theoretically  exist  at  the  points  of 
intersection  of  the  trench  with  the  planar  region  and  which  are  not  approxi¬ 
mated  very  accurately  by  the  mesh  being  used  only,  effect  the  numerical  solu¬ 
tion  in  the  neighborhood  of  the  singularities. 


2.  Steady,  Stably  Stratified  Flew  Over  Topography 


Fig.  ?.  Owens  Valley  Flow 

This  example  illustrates  a  steady  stratified  flow  with  a  Froude  number,  i.e„ 
ratio  of  inertial  to  gravitational  forces,  of  4.5.  For  simplicity  the  simulation  was 
done  specifying  a  constant  ground  temperature  and  an  inlet  temperature  pro¬ 
file  which  increased  linearly  with  height  and  a  specified  inlet  velocity  field  so 
that  a  comparison  could  be  made  with  some  simulations  done  previously  by  the 
group  at  Lawrence  Livermore  National  Laboratory.  The  comparison  was  essen¬ 
tially  exact  and  in  both  cases  such  features  as  mountain  lee  waves  and  up-slope 
winds  are  evident. 

3.  Transient  Neutral  and  Stratified  Flow  About  a  Cylindrical  Object 
* '  This  flow  was  chosen  because  of  the  numer  of  experimental  observations 
available  in  the  literature,  at  least  in  the  neutral  flow  case,  for  comparison.  The 
details  of  the  mesh  and  boundary  conditions  are  displayed  in  Figure  1.  (Note 
that  the  outflow  bounding  conditions  illustrated  are  those  appropriate  to  the 
neutral  flow  case.) 
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&.  Owrall  leau,  atih,  tad  boundary  condition*. 


b.  Douilt  of  aotb  notr  cylinder;  alio  iboon  art  the  nodes  [or  a 
typical  eleaent. 


Fig.  8.  Mesh  and  Boundary  Conditions 

The  numerical  experiment  started  with  the  zero  flow  initially  and  the  velocity 
field  evolves  to  one  at  a  Reynold's  number  of  «  100  (based  on  cylinder  diame¬ 
ter)  with  a  period  shedding  of  a  vortice,  first  from  the  top  and  then  from  the 
bottom  of  the  cylinder.  That  is,  a  von  Karmen  vortex  street  develops  behind  the 
cylinder  and  flows  out  the  computational  domain.  These  features  are  illus¬ 
trated  in  Figures  10-14.  These  predictions  are  in  close  agreement  with  observa¬ 
tions  made  in  physical  experiments. 


Fig.  10.  Details  of  Vortex  Generation 
During  Quarter  Cycle  of  Vortex 
Formation 


Fig.  11.  Time  History  of 
Velocity  At  Node  B 
During  Shedding  Cycle 

The  latter  is  also  illustrated  in  Figures  (12)  -  (13)  where  streaklines  computed 
over  a  numerically  simulated  vortex  shedding  cycle  (Figure  12)  can  be  com¬ 
pared  to  those  experimentally  observed  at  a  Reynolds  number  of  102,  (Figure 
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Fig.  13.  Experimentally  Observed  Streakline 
In  the  comparison  one  can  focus  on  the  "eye  patterns"  as  they  are  swept  down¬ 
stream.  The  downstream  movement  of  these  patterns  away  from  the  centerline 
obvious  in  both  the  numerical  and  physical  experiments  and  overall  features 
compare  very  well.  (The  algorithm  for  computing  the  streaklines  was  developed 
during  the  course  of  our  study  and  has  been  an  aid  in  understanding  some 
features  of  the  flows  as  well  as  a  starting  point  for  tracking  some  neutrally 
buoyant  particles.) 

Finally  to  test  the  new  outflow  boundary  conditions  as  well  as  to  gain  some 
insight  into  a  stratified  flow  about  an  object  the  temperature  of  the  cylinder 
was  raised  above  that  of  the  entering  stream  and  the  top  and  bottom  boun¬ 
daries  were  insulated.  In  our  initial  experiment  the  heating  was  minimal.  The 
overall  flow  structure  didn’t  change  much  except  now  the  vortices  leaving  the 
cylinder  are  more  asymmetric  and  have  a  tendency  to  raise  during  their  down¬ 
stream  course.  Moreover,  the  temperature  field  exhibits  temporal  oscillations 
as  displayed  in  Figure  14. 
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Pig.  14.  Time  History  of  Temperature 
At  Nodes  B,  C  and  D  During 
Shedding  Cycle. 

The  new  outflow  boundary  condition  appeared  to  perform  satisfactory  with 
regard  to  the  velocity  field  calculations  but  some  difficulty  was  encountered  in 
recovering  an  accurate  pressure  field;  the  latter  has  stimulated  some  ongoing 


research  in  this  area. 
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